% %%
% %schematic
clc;clear;
lats=40.5;late=64;lons=135;lone=165;interval=0.0125;
path(path,'/Users/yuannan/Desktop/sea ice in the okhotsk/in japan/m_map1.4/m_map');
c=m_etopo2([lons,lone lats,late]);
[m,n]=size(c);
x4=linspace(lats,late,m);y4=linspace(lons,lone,n);
[x5,y5]=meshgrid(y4,x4);
c=flipud(c);
c(c>0)=nan;
load('lon_int_all.mat');
load('lat_int_all.mat');
c_new=griddata(x5,y5,c,lon_int_all,lat_int_all);
clear lats late lons lone interval

depth=parameter(1);
depth=flipud(depth(1,1:84)')';
lon=parameter(2);lat=parameter(3);
y=size(lon);x=size(lat);y=y(1,2);x=x(1,2);
layer_num=84;interval1=0.1;interval=10;dd=10;
p(1:x,1:y)=10.1325; %reference pressur e dbar

%%
f_id=13;
[fid1,fid2,fid3,fid4,fid5,fid6,fid7,fid8,fid9,tit]=read_data(f_id); 
for layer=84
    fseek(fid1,y*x*4*(layer_num-layer),-1);
    fseek(fid2,y*x*4*(layer_num-layer),-1);
    t=fread(fid1,y*x,'real*4');t1=reshape(t,y,x);t1(t1<-10000)=nan;t_v(:,:,layer)=t1';
    s=fread(fid2,y*x,'real*4');s1=reshape(s,y,x);s1(s1<-10000)=nan;s_v(:,:,layer)=s1';
    rho(:,:,layer)=gsw_rho(s_v(:,:,layer),t_v(:,:,layer),p)-1000; %sa:g/kg ct degc
end
clear t s t1 s1 t_v s_v
fclose(fid1);fclose(fid2);

rho(:,:,1)=rho(:,:,84);
for layer=1:84
    fseek(fid4,y*x*4*(layer-1),-1);
    fseek(fid5,y*x*4*(layer-1),-1);
    u=fread(fid4,y*x,'real*4');u1=reshape(u,y,x);u1(u1<-10000)=nan;u_v(:,:,layer)=u1';
    v=fread(fid5,y*x,'real*4');v1=reshape(v,y,x);v1(v1<-10000)=nan;v_v(:,:,layer)=v1';
end

fclose(fid4);fclose(fid5);
clear u u1 v v1

 %vertical integration and average: method1：先求垂向平均再求meridional平均
range1=1651:1951;
range2=1201:1381;
u_sum_1(1:x,1:y)=0;
v_sum_1(1:x,1:y)=0;
u_ave_1(1:x,1:y)=0;
v_ave_1(1:x,1:y)=0;
num_ver_1(1:x,1:y)=1;
 for ii=range1
     for jj=range2
         for kk=2:layer_num
             if isnan(u_v(ii,jj,kk))
                 continue
             else
                 u_sum_2(ii,jj)=(u_v(ii,jj,kk-1)+u_v(ii,jj,kk))./2.0.*(depth(1,kk)-depth(1,kk-1));
                 v_sum_2(ii,jj)=(v_v(ii,jj,kk-1)+v_v(ii,jj,kk))./2.0.*(depth(1,kk)-depth(1,kk-1));                 
                 u_sum_1(ii,jj)=u_sum_1(ii,jj)+u_sum_2(ii,jj);
                 v_sum_1(ii,jj)=v_sum_1(ii,jj)+v_sum_2(ii,jj);                 
                 num_ver_1(ii,jj)=num_ver_1(ii,jj)+1;
             end
         end
         u_ave_1(ii,jj)=u_sum_1(ii,jj)./depth(1,num_ver_1(ii,jj));
         v_ave_1(ii,jj)=v_sum_1(ii,jj)./depth(1,num_ver_1(ii,jj));
     end
 end
%%
fid_geo=fopen('coco_geo_grid.grd','r');
lon_old=100:0.1:270-0.1;lat_old=-10:0.1:67-0.1;
lon_new=100:1.0/30.0:270-1.0/30.0;lat_new=-10:1.0/30.0:67-1.0/30.0;
[l_old,t_old]=meshgrid(lon_old,lat_old); 
[l_new,t_new]=meshgrid(lon_new,lat_new);
c2_old=fread(fid_geo,1700*770,'real*4');
c2_old=reshape(c2_old,1700,770);c2_old=c2_old';
c2=interp2(l_old,t_old,c2_old,l_new,t_new);
c2(c2<-100000)=nan;

%%
%wind stress
srcDir1=uigetdir('/Users/yuannan/Desktop/sea ice in the okhotsk/ERA-japan/40years_data/original_data_program/monthly_mean_39_years_ECMWF');
cd(srcDir1);
lons=0;lone=359.25;lats=90;late=-90;
iter=0.75;lonc=lons:iter:lone;latc=lats:-1*iter:late;
ix=size(lonc);iy=size(latc);ix=ix(1,2);iy=iy(1,2);
[LON_C, LAT_C]=meshgrid(lonc,latc);
fid1=fopen('TAUX_monthly_39_real.grd','r');fid2=fopen('TAUY_monthly_39_real.grd','r');
taux=fread(fid1,ix*iy,'real*4');
tauy=fread(fid2,ix*iy,'real*4');
taux=reshape(taux,ix,iy);
tauy=reshape(tauy,ix,iy);
taux_new=interp2(LON_C,LAT_C,taux',l_new,t_new);
tauy_new=interp2(LON_C,LAT_C,tauy',l_new,t_new);

%%
% range_c1=31:871;range_c2=301:901;%[41 55];[140 150]
% % sf=surf(x5(range_c1,range_c2),y5(range_c1,range_c2),c(range_c1,range_c2));
% [fig_3dx,fig_3dy]=meshgrid(1:800,1:800);
% sf=surf(fig_3dx,fig_3dy,c_new');colormap(othercolor('Accent5'));
% set(sf,'edgecolor','none','facecolor','interp');
% view(75,69);caxis([-2500 0]);
% hold on;contour3(fig_3dx(1:420,1:350),fig_3dy(1:420,1:350),c_new(1:350,1:420)',[-50 -100 -200 -500 -1000],'linewidth',2,'color','k');
% set(gca,'zlim',[-1800 0],'xlim',[100 412],'ylim',[200 400],'xtick',[],'ytick',[],'ztick',[],'xcolor','none','ycolor','none','zcolor','none')
% hold on;[c21,h21]=contour(fig_3dx,fig_3dy,lon_int_all',[140 141 142 143 144 145 146 147 148 149 150 151 152 153],'linewidth',2,'color',[.7 .7 .7],'linestyle',':');
% hold on;[c11,h11]=contour(fig_3dx,fig_3dy,lat_int_all',[44 45 46 47 48 49 50 51 52 53 54 55 56],'linewidth',2,'color',[.7 .7 .7],'linestyle',':');
% set(gcf,'unit','centimeters','position',[0,0,80,80]);
% print('topo_new1','-dtiff','-r600');
%%
cd /Users/yuannan/Desktop/paper_figures
pcolor(x5,y5,-1*c);shading interp;colormap(flipud(cmocean('deep')));
hold on;contour(l_new,t_new,c2,[0 0],'linewidth',3,'color','k');
hold on;[ct,ht]=contour(l_new(1621:end,:),t_new(1621:end,:),c2(1621:end,:),[50 50],'linewidth',3,'color',[255 102 0]/255);
hold on;[ct,ht]=contour(l_new(1621:end,:),t_new(1621:end,:),c2(1621:end,:),[100 100],'linewidth',3,'color',[245 222 179]/255);
hold on;[ct,ht]=contour(l_new(1621:end,:),t_new(1621:end,:),c2(1621:end,:),[200 200],'linewidth',3,'color',[255 215 0]/255);
hold on;[ct,ht]=contour(l_new(1621:end,:),t_new(1621:end,:),c2(1621:end,:),[500 500],'linewidth',3,'color',[255 165 0]/255);
hold on;[ct,ht]=contour(l_new(1621:end,:),t_new(1621:end,:),c2(1621:end,:),[1000 1000],'linewidth',3,'color',[255 250 205]/255);

% hold on;[ct,ht]=contour(l_new(1621:end,:),t_new(1621:end,:),c2(1621:end,:),[50 50],'linewidth',3,'color',[100 25 25]/255);
% hold on;[ct,ht]=contour(l_new(1621:end,:),t_new(1621:end,:),c2(1621:end,:),[100 100],'linewidth',3,'color',[100 20 5]/255);
% hold on;[ct,ht]=contour(l_new(1621:end,:),t_new(1621:end,:),c2(1621:end,:),[200 200],'linewidth',3,'color',[80 10 0]/255);
% hold on;[ct,ht]=contour(l_new(1621:end,:),t_new(1621:end,:),c2(1621:end,:),[500 500],'linewidth',3,'color',[95 65 30]/255);
% hold on;[ct,ht]=contour(l_new(1621:end,:),t_new(1621:end,:),c2(1621:end,:),[1000 1000],'linewidth',3,'color',[90 60 0]/255);
% 
%clabel(ct,ht,'fontsize',10,'labelspacing',600,'FontName','Arial Rounded MT Bold','color',[165 42 42]/255);
caxis([-3000 7000]);
colb=colorbar;set(colb,'limits',[0 7000],'Ticks',[0 2000 4000 6000],'TickLabels',{'0','2000','4000','6000'},'fontsize',15);
title(colb,'depth(m)','FontName','Arial Rounded MT Bold','fontsize',15);
hold on;line([141 141 145 145 141],[55 52 52 55 55],'linewidth',4,'color','k');
hold on;line([141 141 145 145 141],[55 52 52 55 55],'linewidth',4,'color',[255 105 180]/255,'linestyle','--');
set(gca,'linewidth',4,'Yaxislocation','left','tickdir','out','FontName','Arial Rounded MT Bold','fontsize',25);
set(gca,'ylim',[41 56],'ytick',[41 44 47 50 53 56],'xlim',[138 150],'xtick',[138 141 144 147 150]);
set(gca,'yticklabel',{'41^oN','44^oN','47^oN','50^oN','53^oN','56^oN'},'xticklabel',{'138^oE','141^oE','144^oE','147^oE','150^oE'});
%text(139,53.3,'Amur River','FontName','Arial Rounded MT Bold','fontsize',15,'rotation',-15);
size1=size(y5);
[xx5,yy5]=meshgrid(1:800,1:800);
xx5_raw=griddata(lon_int_all,lat_int_all,xx5,x5,y5);
yy5_raw=griddata(lon_int_all,lat_int_all,yy5,x5,y5);
hold on;contour(x5,y5,yy5_raw,[50 100 150 200 250 300 350 400 450],'linewidth',3,'color',[211 211 211]/255);
hold on;contour(x5,y5,xx5_raw,[150 180 210 240 270 300 330 360 390 420 450],'linewidth',3,'color',[211 211 211]/255);
hold on;
quiver(l_new(1:22:end,1:22:end),t_new(1:22:end,1:22:end),taux_new(1:22:end,1:22:end),tauy_new(1:22:end,1:22:end),2,'color',[0,0,139]/255,'linewidth',1.5,'maxheadsize',5);hold on;
hold on;
% quiver(147,42.25,0.2,0,2,'color',[0,0,139]/255,'linewidth',1.5,'maxheadsize',10);
% text(147.5,42.25,'\bf{\tau}(0.2 Pa)','FontName','Arial Rounded MT Bold','fontsize',15,'color',[0,0,139]/255);
quiver(138.5,52,0,0.02,2,'color',[0,0,139]/255,'linewidth',1.5,'maxheadsize',10);
text(139,56.5,'\bf{\tau}(0.2 Pa)','FontName','Arial Rounded MT Bold','fontsize',15,'color',[0,0,139]/255);

% text(142.5,48.7,'Sakhalin  Island','FontName','Arial Rounded MT Bold','fontsize',24,'rotation',90,'color',[255 244 230]/255);
% text(141.6,43.4,{'Hokkaido','  Island'},'FontName','Arial Rounded MT Bold','fontsize',24,'color',[255 244 230]/255);
text(142.5,48.7,'Sakhalin  Island','FontName','Arial Rounded MT Bold','fontsize',24,'rotation',90,'color','k');
text(141.6,43.4,{'Hokkaido','  Island'},'FontName','Arial Rounded MT Bold','fontsize',24,'color','k');

set(gcf,'unit','centimeters','position',[0,0,23,40]);
%print('topo_3','-dtiff','-r600');
%%
cd /Users/yuannan/Desktop/paper_figures
figure(2)
range1=1651:1951;
range2=1201:1381;
range1_1=1865:4:1951;
range2_1=1200:4:1381;
m_proj('miller','lon',[141 145],'lat',[52 55]);
m_pcolor(l_new(range1,range2),t_new(range1,range2),rho(range1,range2));shading interp;colormap(cmocean('dense'));
caxis([20 27]);
hold on;m_contour(l_new,t_new,c2,[0 0],'linewidth',3,'color','k');
hold on;[ct,ht]=m_contour(l_new(1621:end,:),t_new(1621:end,:),c2(1621:end,:),[50 50],'linewidth',3,'color',[255 102 0]/255,'showtext','on','labelspacing',350);
clabel(ct,ht,'FontName','Arial Rounded MT Bold','fontsize',15,'color',[255 102 0]/255);
hold on;[ct,ht]=m_contour(l_new(1621:end,:),t_new(1621:end,:),c2(1621:end,:),[100 100],'linewidth',3,'color',[245 222 179]/255,'showtext','on','labelspacing',350);
clabel(ct,ht,'FontName','Arial Rounded MT Bold','fontsize',15,'color',[245 222 179]/255);
hold on;[ct,ht]=m_contour(l_new(1621:end,:),t_new(1621:end,:),c2(1621:end,:),[200 200],'linewidth',3,'color',[255 215 0]/255,'showtext','on','labelspacing',350);
clabel(ct,ht,'FontName','Arial Rounded MT Bold','fontsize',15,'color',[255 215 0]/255);
hold on;[ct,ht]=m_contour(l_new(1621:end,:),t_new(1621:end,:),c2(1621:end,:),[500 500],'linewidth',3,'color',[255 165 0]/255,'showtext','on','labelspacing',350);
clabel(ct,ht,'FontName','Arial Rounded MT Bold','fontsize',15,'color',[255 165 0]/255);
hold on;[ct,ht]=m_contour(l_new(1621:end,:),t_new(1621:end,:),c2(1621:end,:),[1000 1000],'linewidth',3,'color',[255 250 205]/255,'showtext','on','labelspacing',350);
clabel(ct,ht,'FontName','Arial Rounded MT Bold','fontsize',15,'color',[255 250 205]/255);
hold on;[ct,ht]=m_contour(l_new(1621:end,:),t_new(1621:end,:),rho(1621:end,:,1),[15 17 19 21 23 25 26 27],'linewidth',3,'color','k','showtext','on','labelspacing',280);
clabel(ct,ht,'FontName','Arial Rounded MT Bold','fontsize',15);
hold on;
m_vec(70,l_new(range1_1,range2_1),t_new(range1_1,range2_1),u_ave_1(range1_1,range2_1),v_ave_1(range1_1,range2_1),'c','shaftwidth',1,'headlength',5);
m_grid('tickdir','out','linewi',4,'fontsize',25,'ytick',[52 53 54 55]);
m_gshhs_h('patch',[.7 .7 .7],'linewidth',4);
hold on;
m_vec(70,142,53.1,10,0,'r','shaftwidth',1,'headlength',5);
hold on;
line1(1:10)=53.5;hold on;m_line(143.1:0.1:144,line1,'linewidth',5,'color',[206 251 109]/255);
hold on;
m_text(142.2,53.1,'10cm/s','FontName','Arial Rounded MT Bold','fontsize',13,'color','r');
m_text(142,52.9,'barotropic current','FontName','Arial Rounded MT Bold','fontsize',13,'color','r');
set(gca,'linewidth',4,'ytick',[52 53 54 55],'xtick',[141 142 143 144 145],'FontName','Arial Rounded MT Bold','fontsize',18);
%axis([141 145 52 55]);
col=colorbar;
title(col,'\bf{\rho}(kg/m^3-1000)','FontName','Arial Rounded MT Bold','color',[255 250 205]/255,'fontsize',14);
set(col,'position',[0.8 0.55 0.05 0.3],'color',[255 250 205]/255,'fontsize',14);
% set(gca,'linewidth',4,'Yaxislocation','left','tickdir','out','FontName','Arial Rounded MT Bold','fontsize',25);
% set(gca,'ylim',[52 55],'xlim',[141 145]);
% set(gca,'yticklabel',{'52^oN','53^oN','54^oN','55^oN'},'xticklabel',{'141^oE','142^oE','143^oE','144^oE','145^oE'});
set(gcf,'unit','centimeters','position',[0,0,20,21]);
print('GRL1_section2_new','-dtiff','-r600');

